home *** CD-ROM | disk | FTP | other *** search
/ Digital Information Mana…ntial Guide to Multimedia / Digital Information Management - An Essential Guide to Multimedia.iso / Audacity / Nyquist / dspprims.lsp < prev    next >
Lisp/Scheme  |  2006-05-16  |  17KB  |  521 lines
  1. ;; dspprims.lsp -- interface to dsp primitives
  2.  
  3. ;; ARESON - notch filter
  4. ;; 
  5. (defun areson (s c b &optional (n 0))
  6.   (multichan-expand #'nyq:areson s c b n))
  7.  
  8. (setf areson-implementations
  9.       (vector #'snd-areson #'snd-aresonvc #'snd-aresoncv #'snd-aresonvv))
  10.  
  11. ;; NYQ:ARESON - notch filter, single channel
  12. ;;
  13. (defun nyq:areson (signal center bandwidth normalize)
  14.   (select-implementation-1-2 areson-implementations 
  15.    signal center bandwidth normalize))
  16.  
  17.  
  18. ;; hp - highpass filter
  19. ;; 
  20. (defun hp (s c)
  21.   (multichan-expand #'nyq:hp s c))
  22.  
  23. (setf hp-implementations
  24.       (vector #'snd-atone #'snd-atonev))
  25.  
  26. ;; NYQ:hp - highpass filter, single channel
  27. ;;
  28. (defun nyq:hp (s c)
  29.   (select-implementation-1-1 hp-implementations s c))
  30.  
  31.  
  32. ;; comb-delay-from-hz -- compute the delay argument
  33. ;;
  34. (defun comb-delay-from-hz (hz caller)
  35.   (recip hz))
  36.  
  37. ;; comb-feedback-from-decay -- compute the feedback argument
  38. ;;
  39. (defun comb-feedback (decay delay)
  40.   (s-exp (mult -6.9087 delay (recip decay))))
  41.  
  42. ;; COMB - comb filter
  43. ;; 
  44. ;; this is just a feedback-delay with different arguments
  45. ;;
  46. (defun comb (snd decay hz)
  47.   (multichan-expand #'nyq:comb snd decay hz))
  48.  
  49. (defun nyq:comb (snd decay hz)
  50.   (let (delay feedback len d)
  51.     ; convert decay to feedback, iterate over array if necessary
  52.     (setf delay (comb-delay-from-hz hz "comb"))
  53.     (setf feedback (comb-feedback decay delay))
  54.     (nyq:feedback-delay snd delay feedback)))
  55.  
  56. ;; ALPASS - all-pass filter
  57. ;; 
  58. (defun alpass (snd decay hz &optional min-hz)
  59.   (multichan-expand #'nyq:alpass snd decay hz min-hz))
  60.   
  61.  
  62.  
  63. (defun nyq:alpass (snd decay hz min-hz)
  64.   (let (delay feedback len d)
  65.     ; convert decay to feedback, iterate over array if necessary
  66.     (setf delay (comb-delay-from-hz hz "alpass"))
  67.     (setf feedback (comb-feedback decay delay))
  68.     (nyq:alpass1 snd delay feedback min-hz)))
  69.  
  70.  
  71. ;; CONST -- a constant at control-srate
  72. ;;
  73. (defun const (value &optional (dur 1.0))
  74.   (let ((d (get-duration dur)))
  75.     (snd-const value *rslt* *CONTROL-SRATE* d)))
  76.  
  77.  
  78. ;; CONVOLVE - slow convolution
  79. ;; 
  80. (defun convolve (s r)
  81.   (multichan-expand #'snd-convolve s r))
  82.  
  83.  
  84. ;; FEEDBACK-DELAY -- (delay is quantized to sample period)
  85. ;;
  86. (defun feedback-delay (snd delay feedback)
  87.   (multichan-expand #'nyq:feedback-delay snd delay feedback))
  88.   
  89.  
  90. ;; SND-DELAY-ERROR -- report type error
  91. ;;
  92. (defun snd-delay-error (snd delay feedback)
  93.   (error "feedback-delay with variable delay is not implemented"))
  94.  
  95.  
  96. ;; NYQ::DELAYCV -- coerce sample rates and call snd-delaycv
  97. ;;
  98. (defun nyq:delaycv (the-snd delay feedback)
  99.   (display "delaycv" the-snd delay feedback)
  100.   (let ((the-snd-srate (snd-srate the-snd))
  101.         (feedback-srate (snd-srate feedback)))
  102.     (cond ((> the-snd-srate feedback-srate)
  103.            (setf feedback (snd-up the-snd-srate feedback)))
  104.           ((< the-snd-srate feedback-srate)
  105.            (format t "Warning: down-sampling feedback in feedback-delay/comb~%")
  106.            (setf feedback (snd-down the-snd-srate feedback))))
  107.     (snd-delaycv the-snd delay feedback)))
  108.  
  109. (setf feedback-delay-implementations
  110.       (vector #'snd-delay #'snd-delay-error #'nyq:delaycv #'snd-delay-error))
  111.  
  112.  
  113. ;; NYQ:FEEDBACK-DELAY -- single channel delay
  114. ;;
  115. (defun nyq:feedback-delay (snd delay feedback)
  116.   (select-implementation-1-2 feedback-delay-implementations 
  117.                              snd delay feedback))
  118.  
  119.  
  120. ;; SND-ALPASS-ERROR -- report type error
  121. ;;
  122. (defun snd-alpass-error (snd delay feedback)
  123.   (error "alpass with constant decay and variable hz is not implemented"))
  124.  
  125.  
  126. (if (not (fboundp 'snd-alpasscv))
  127.     (defun snd-alpasscv (snd delay feedback min-hz)
  128.       (error "snd-alpasscv (ALPASS with variable decay) is not implemented")))
  129. (if (not (fboundp 'snd-alpassvv))
  130.     (defun snd-alpassvv (snd delay feedback min-hz)
  131.       (error "snd-alpassvv (ALPASS with variable decay and feedback) is not implemented")))
  132.       
  133. (defun snd-alpass-4 (snd delay feedback min-hz)
  134.     (snd-alpass snd delay feedback))
  135.     
  136.  
  137. (defun snd-alpasscv-4 (the-snd delay feedback min-hz)
  138.     (display "snd-alpasscv-4" (snd-srate the-snd) (snd-srate feedback))
  139.     (let ((the-snd-srate (snd-srate the-snd))
  140.           (feedback-srate (snd-srate feedback)))
  141.       (cond ((> the-snd-srate feedback-srate)
  142.              (setf feedback (snd-up the-snd-srate feedback)))
  143.             ((< the-snd-srate feedback-srate)
  144.              (format t "Warning: down-sampling feedback in alpass~%")
  145.              (setf feedback (snd-down the-snd-srate feedback))))
  146.       (display "snd-alpasscv-4 after cond" (snd-srate the-snd) (snd-srate feedback))
  147.       (snd-alpasscv the-snd delay feedback)))
  148.  
  149.     
  150. (defun snd-alpassvv-4 (the-snd delay feedback min-hz)
  151.     ;(display "snd-alpassvv-4" (snd-srate the-snd) (snd-srate feedback))
  152.     (let ((the-snd-srate (snd-srate the-snd))
  153.           (delay-srate (snd-srate delay))
  154.           (feedback-srate (snd-srate feedback))
  155.           max-delay)
  156.       (cond ((or (not (numberp min-hz))
  157.                  (<= min-hz 0))
  158.              (error "alpass needs numeric (>0) 4th parameter (min-hz) when delay is variable")))
  159.       (setf max-delay (/ 1.0 min-hz))
  160.       ; make sure delay is between 0 and max-delay
  161.       ; use clip function, which is symetric, with an offset
  162.       (setf delay (snd-offset (clip (snd-offset delay (* max-delay 0.5))
  163.                                     max-delay)
  164.                               (* max-delay -0.5)))
  165.       ; now delay is between 0 and max-delay, so we won't crash nyquist when
  166.       ; we call snd-alpassvv, which doesn't test for out-of-range data
  167.       (cond ((> the-snd-srate feedback-srate)
  168.              (setf feedback (snd-up the-snd-srate feedback)))
  169.             ((< the-snd-srate feedback-srate)
  170.              (format t "Warning: down-sampling feedback in alpass~%")
  171.              (setf feedback (snd-down the-snd-srate feedback))))
  172.       (cond ((> the-snd-srate delay-srate)
  173.              (setf delay (snd-up the-snd-srate delay)))
  174.             ((< the-snd-srate delay-srate)
  175.              (format t "Warning: down-sampling delay in alpass~%")
  176.              (setf delay (snd-down the-snd-srate delay))))
  177.       ;(display "snd-alpassvv-4 after cond" (snd-srate the-snd) (snd-srate feedback))
  178.       (snd-alpassvv the-snd delay feedback max-delay)))
  179.  
  180. (setf alpass-implementations
  181.       (vector #'snd-alpass-4 #'snd-alpass-error
  182.               #'snd-alpasscv-4 #'snd-alpassvv-4))
  183.  
  184.  
  185.  
  186. ;; NYQ:ALPASS1 -- single channel alpass
  187. ;;
  188. (defun nyq:alpass1 (snd delay feedback min-hz)
  189.   (select-implementation-1-2 alpass-implementations
  190.                              snd delay feedback min-hz))
  191.  
  192. ;; S-EXP -- exponentiate a sound
  193. ;;
  194. (defun s-exp (s) (multichan-expand #'nyq:exp s))
  195.  
  196.  
  197. ;; NYQ:EXP -- exponentiate number or sound
  198. ;;
  199. (defun nyq:exp (s) (if (soundp s) (snd-exp s) (exp s)))
  200.  
  201. ;; S-ABS -- absolute value of a sound
  202. ;;
  203. (defun s-abs (s) (multichan-expand #'nyq:abs s))
  204.  
  205. ;; NYQ:ABS -- absolute value of number or sound
  206. ;;
  207. (defun nyq:abs (s) (if (soundp s) (snd-abs s) (abs s)))
  208.  
  209. ;; S-SQRT -- square root of a sound
  210. ;;
  211. (defun s-sqrt (s) (multichan-expand #'nyq:sqrt s))
  212.  
  213. ;; NYQ:SQRT -- square root of a number or sound
  214. ;;
  215. (defun nyq:sqrt (s) (if (soundp s) (snd-sqrt s) (sqrt s)))
  216.  
  217.  
  218. ;; INTEGRATE -- integration
  219. ;;
  220. (defun integrate (s) (multichan-expand #'snd-integrate s))
  221.  
  222.  
  223. ;; S-LOG -- natural log of a sound
  224. ;;
  225. (defun s-log (s) (multichan-expand #'nyq:log s))
  226.  
  227.  
  228. ;; NYQ:LOG -- log of a number or sound
  229. ;;
  230. (defun nyq:log (s) (if (soundp s) (snd-log s) (log s)))
  231.  
  232.  
  233. ;; NOISE -- white noise
  234. ;;
  235. (defun noise (&optional (dur 1.0))
  236.   (let ((d (get-duration dur)))
  237.     (snd-white *rslt* *SOUND-SRATE* d)))
  238.  
  239.  
  240. (defun noise-gate (snd &optional (lookahead 0.5) (risetime 0.02) (falltime 0.5)
  241.                                                  (floor 0.01) (threshold 0.01))
  242.   (let ((rms (lp (mult snd snd) (/ *control-srate* 10.0))))
  243.     (setf threshold (* threshold threshold))
  244.     (mult snd (gate rms lookahead risetime falltime floor threshold))))
  245.  
  246.  
  247. ;; QUANTIZE -- quantize a sound
  248. ;;
  249. (defun quantize (s f) (multichan-expand #'snd-quantize s f))
  250.  
  251.  
  252. ;; RECIP -- reciprocal of a sound
  253. ;;
  254. (defun recip (s) (multichan-expand #'nyq:recip s))
  255.  
  256.  
  257. ;; NYQ:RECIP -- reciprocal of a number or sound
  258. ;;
  259. (defun nyq:recip (s) (if (soundp s) (snd-recip s) (/ (float s))))
  260.  
  261. ;; RMS -- compute the RMS of a sound
  262. ;;
  263. (defun rms (s &optional (rate 100.0) window-size)
  264.   (let (rslt step-size)
  265.     (setf step-size (round (/ (snd-srate s) rate)))
  266.     (cond ((null window-size)
  267.                (setf window-size step-size)))
  268.     (setf s (prod s s))
  269.     (setf result (snd-avg s window-size step-size OP-AVERAGE))
  270.         ;; compute square root of average
  271.         (s-exp (scale 0.5 (s-log result)))))
  272.  
  273.  
  274. ;; RESON - bandpass filter
  275. ;; 
  276. (defun reson (s c b &optional (n 0))
  277.   (multichan-expand #'nyq:reson s c b n))
  278.  
  279. (setf reson-implementations
  280.       (vector #'snd-reson #'snd-resonvc #'snd-resoncv #'snd-resonvv))
  281.  
  282. ;; NYQ:RESON - bandpass filter, single channel
  283. ;;
  284. (defun nyq:reson (signal center bandwidth normalize)
  285.   (select-implementation-1-2 reson-implementations 
  286.    signal center bandwidth normalize))
  287.  
  288.  
  289. ;; SHAPE -- waveshaper
  290. ;;
  291. (defun shape (snd shape origin)
  292.   (multichan-expand #'snd-shape snd shape origin))
  293.  
  294.  
  295. ;; SLOPE -- calculate the first derivative of a signal
  296. ;;
  297. (defun slope (s) (multichan-expand #'nyq:slope s))
  298.  
  299.  
  300. ;; NYQ:SLOPE -- first derivative of single channel
  301. ;;
  302. (defun nyq:slope (s)
  303.   (let* ((sr (snd-srate s))
  304.          (sr-inverse (/ sr)))
  305.     (snd-xform (snd-slope s) sr (- sr-inverse) 0.0 MAX-STOP-TIME 1.0)))
  306.  
  307.  
  308. ;; lp - lowpass filter
  309. ;; 
  310. (defun lp (s c)
  311.   (multichan-expand #'nyq:lp s c))
  312.  
  313. (setf lp-implementations
  314.       (vector #'snd-tone #'snd-tonev))
  315.  
  316. ;; NYQ:lp - lowpass filter, single channel
  317. ;;
  318. (defun nyq:lp (s c)
  319.   (select-implementation-1-1 lp-implementations s c))
  320.  
  321.  
  322.  
  323. ;;; fixed-parameter filters based on snd-biquad
  324.  
  325. (setf Pi 3.14159265358979)
  326.  
  327. (defun square (x) (* x x))
  328. (defun sinh (x) (* 0.5 (- (exp x) (exp (- x)))))
  329.  
  330.  
  331. ; remember that snd-biquad uses the opposite sign convention for a_i's 
  332. ; than Matlab does.
  333.  
  334. ; convenient biquad: normalize a0, and use zero initial conditions.
  335. (defun biquad (x b0 b1 b2 a0 a1 a2)
  336.   (let ((a0r (/ 1.0 a0)))
  337.     (snd-biquad x (* a0r b0) (* a0r b1) (* a0r b2) 
  338.                              (* a0r a1) (* a0r a2) 0 0)))
  339.  
  340. ; biquad with Matlab sign conventions for a_i's.
  341. (defun biquad-m (x b0 b1 b2 a0 a1 a2)
  342.   (biquad x b0 b1 b2 a0 (- a1) (- a2)))
  343.  
  344. ; two-pole lowpass
  345. (defun lowpass2 (x hz &optional (q 0.7071))
  346.   (let* ((w (* 2.0 Pi (/ hz (snd-srate x))))
  347.          (cw (cos w))
  348.          (sw (sin w))
  349.          (alpha (* sw (sinh (/ 0.5 q))))
  350.          (a0 (+ 1.0 alpha))
  351.          (a1 (* -2.0 cw))
  352.          (a2 (- 1.0 alpha))
  353.          (b1 (- 1.0 cw))
  354.          (b0 (* 0.5 b1))
  355.          (b2 b0))
  356.     (biquad-m x b0 b1 b2 a0 a1 a2)))
  357.  
  358. ; two-pole highpass
  359. (defun highpass2 (x hz &optional (q 0.7071))
  360.   (let* ((w (* 2.0 Pi (/ hz (snd-srate x))))
  361.          (cw (cos w))
  362.          (sw (sin w))
  363.          (alpha (* sw (sinh (/ 0.5 q))))
  364.          (a0 (+ 1.0 alpha))
  365.          (a1 (* -2.0 cw))
  366.          (a2 (- 1.0 alpha))
  367.          (b1 (- -1.0 cw))
  368.          (b0 (* -0.5 b1))
  369.          (b2 b0))
  370.     (biquad-m x b0 b1 b2 a0 a1 a2)))
  371.  
  372. ; two-pole bandpass.  max gain is unity.
  373. (defun bandpass2 (x hz q)
  374.   (let* ((w (* 2.0 Pi (/ hz (snd-srate x))))
  375.          (cw (cos w))
  376.          (sw (sin w))
  377.          (alpha (* sw (sinh (/ 0.5 q))))
  378.          (a0 (+ 1.0 alpha))
  379.          (a1 (* -2.0 cw))
  380.          (a2 (- 1.0 alpha))
  381.          (b0 alpha)
  382.          (b1 0.0)
  383.          (b2 (- alpha)))
  384.     (biquad-m x b0 b1 b2 a0 a1 a2)))
  385.  
  386. ; two-pole notch.
  387. (defun notch2 (x hz q)
  388.   (let* ((w (* 2.0 Pi (/ hz (snd-srate x))))
  389.          (cw (cos w))
  390.          (sw (sin w))
  391.          (alpha (* sw (sinh (/ 0.5 q))))
  392.          (a0 (+ 1.0 alpha))
  393.          (a1 (* -2.0 cw))
  394.          (a2 (- 1.0 alpha))
  395.          (b0 1.0)
  396.          (b1 (* -2.0 cw))
  397.          (b2 1.0))
  398.     (biquad-m x b0 b1 b2 a0 a1 a2)))
  399.  
  400.  
  401. ; two-pole allpass.
  402. (defun allpass2 (x hz q)
  403.   (let* ((w (* 2.0 Pi (/ hz (snd-srate x))))
  404.          (cw (cos w))
  405.          (sw (sin w))
  406.          (k (exp (* -0.5 w (/ 1.0 q))))
  407.          (a0 1.0)
  408.          (a1 (* -2.0 cw k))
  409.          (a2 (* k k))
  410.          (b0 a2)
  411.          (b1 a1)
  412.          (b2 1.0))
  413.     (biquad-m x b0 b1 b2 a0 a1 a2)))
  414.  
  415.  
  416. ; bass shelving EQ.  gain in dB; Fc is halfway point.
  417. ; response becomes peaky at slope > 1.
  418. (defun eq-lowshelf (x hz gain &optional (slope 1.0))
  419.   (let* ((w (* 2.0 Pi (/ hz (snd-srate x))))
  420.          (sw (sin w))
  421.          (cw (cos w))
  422.          (A (expt 10.0 (/ gain (* 2.0 20.0))))
  423.          (b (sqrt (- (/ (+ 1.0 (square A)) slope) (square (- A 1.0)))))
  424.          (apc (* cw (+ A 1.0)))
  425.          (amc (* cw (- A 1.0)))
  426.          (bs (* b sw))
  427.  
  428.          (b0 (*      A (+ A  1.0 (- amc)    bs  )))
  429.          (b1 (*  2.0 A (+ A -1.0 (- apc)        )))
  430.          (b2 (*      A (+ A  1.0 (- amc) (- bs) )))
  431.          (a0           (+ A  1.0    amc     bs  ))
  432.          (a1 (* -2.0   (+ A -1.0    apc         )))
  433.          (a2           (+ A  1.0    amc  (- bs) )))
  434.     (biquad-m x b0 b1 b2 a0 a1 a2)))
  435.  
  436.  
  437. ; treble shelving EQ.  gain in dB; Fc is halfway point.
  438. ; response becomes peaky at slope > 1.
  439. (defun eq-highshelf (x hz gain &optional (slope 1.0))
  440.   (let* ((w (* 2.0 Pi (/ hz (snd-srate x))))
  441.          (sw (sin w))
  442.          (cw (cos w))
  443.          (A (expt 10.0 (/ gain (* 2.0 20.0))))
  444.          (b (sqrt (- (/ (+ 1.0 (square A)) slope) (square (- A 1.0)))))
  445.          (apc (* cw (+ A 1.0)))
  446.          (amc (* cw (- A 1.0)))
  447.          (bs (* b sw))
  448.  
  449.          (b0 (*      A (+ A  1.0    amc     bs  )))
  450.          (b1 (* -2.0 A (+ A -1.0    apc         )))
  451.          (b2 (*      A (+ A  1.0    amc  (- bs) )))
  452.          (a0           (+ A  1.0 (- amc)    bs  ))
  453.          (a1 (*  2.0   (+ A -1.0 (- apc)        )))
  454.          (a2           (+ A  1.0 (- amc) (- bs) )))
  455.     (biquad-m x b0 b1 b2 a0 a1 a2)))
  456.     
  457. (defun nyq:eq-band (x hz gain width)
  458.   (cond ((and (numberp hz) (numberp gain) (numberp width))
  459.          (eq-band-ccc x hz gain width))
  460.         ((and (soundp hz) (soundp gain) (soundp width))
  461.          (snd-eqbandvvv x hz (db-to-linear gain) width))
  462.         (t
  463.          (error "eq-band hz, gain, and width must be all numbers or all sounds"))))
  464.  
  465. ; midrange EQ.  gain in dB, width in octaves (half-gain width).
  466. (defun eq-band (x hz gain width)
  467.   (multichan-expand #'nyq:eq-band x hz gain width))
  468.   
  469.   
  470. (defun eq-band-ccc (x hz gain width)
  471.   (let* ((w (* 2.0 Pi (/ hz (snd-srate x))))
  472.          (sw (sin w))
  473.          (cw (cos w))
  474.          (J (sqrt (expt 10.0 (/ gain 20.0))))
  475.          ;(dummy (display "eq-band-ccc" gain J))
  476.          (g (* sw (sinh (* 0.5 (log 2.0) width (/ w sw)))))
  477.          ;(dummy2 (display "eq-band-ccc" width w sw g))
  478.          (b0 (+ 1.0 (* g J)))
  479.          (b1 (* -2.0 cw))
  480.          (b2 (- 1.0 (* g J)))
  481.          (a0 (+ 1.0 (/ g J)))
  482.          (a1 (- b1))
  483.          (a2 (- (/ g J) 1.0)))
  484.     (biquad x b0 b1 b2 a0 a1 a2)))
  485.  
  486. ; see failed attempt in eub-reject.lsp to do these with higher-order fns:
  487.  
  488. ; four-pole Butterworth lowpass
  489. (defun lowpass4 (x hz)
  490.   (lowpass2 (lowpass2 x hz 0.60492333) hz 1.33722126))
  491.  
  492. ; six-pole Butterworth lowpass
  493. (defun lowpass6 (x hz)
  494.   (lowpass2 (lowpass2 (lowpass2 x hz 0.58338080) 
  495.                                   hz 0.75932572) 
  496.                                   hz 1.95302407))
  497.  
  498. ; eight-pole Butterworth lowpass
  499. (defun lowpass8 (x hz)
  500.   (lowpass2 (lowpass2 (lowpass2 (lowpass2 x hz 0.57622191)
  501.                                             hz 0.66045510) 
  502.                                             hz 0.94276399)
  503.                                             hz 2.57900101))
  504.  
  505. ; four-pole Butterworth highpass
  506. (defun highpass4 (x hz)
  507.   (highpass2 (highpass2 x hz 0.60492333) hz 1.33722126))
  508.  
  509. ; six-pole Butterworth highpass
  510. (defun highpass6 (x hz)
  511.   (highpass2 (highpass2 (highpass2 x hz 0.58338080) 
  512.                                      hz 0.75932572) 
  513.                                      hz 1.95302407))
  514.  
  515. ; eight-pole Butterworth highpass
  516. (defun highpass8 (x hz)
  517.   (highpass2 (highpass2 (highpass2 (highpass2 x hz 0.57622191)
  518.                                                 hz 0.66045510) 
  519.                                                 hz 0.94276399)
  520.                                                 hz 2.57900101))
  521.